Quantum fluctuation theorem for heat exchange in the strong coupling regime 
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We study quantum heat exchange in a multi-state impurity coupled to two thermal reservoirs. 
Allowing for strong system-bath interactions, we show that a steady-state heat exchange fluctuation 
theorem holds, though the dynamical processes nonlinearly involve the two reservoirs. We accom- 
plish a closed expression for the cumulant generating function, and use it obtain the heat current 
and its cumulants in a nonlinear thermal junction, the two-bath spin-boson model. 
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Exact fluctuation relations for nonequilibrium classical 
systems have been recently discovered and exemplified, 
dealing with work and entropy fluctuations in various 
(open, closed, driven) systems [1]. In particular, the fluc- 
tuation theorem (FT) for entropy production quantifies 
the probability of negative entropy generation, measur- 
ing "second law violation" [2, 3]. Both transient and 
steady-state fluctuation theorems (SSFT) have been de- 
rived, where the latter measures entropy production in 
nonequilibrium steady-state systems over a long interval. 
In the context of heat exchange between two equilibrium 
reservoirs, v = L,R, the SSFT can be roughly stated as 
ln[P t (+u))/P t (-u))] = A/3w [4, 5]. Here P t (u) denotes the 
probability distribution of the net heat transfer w, from 
L to R, over the (long) interval t, with A/3 — T^ 1 - T^ 1 
as the difference between the inverse temperatures of the 
reservoirs. A related quantity is the cumulant generat- 
ing function (CGF), providing general relations between 
transport coefficients under the FT symmetry [6, 7]. 

Extending the work and heat FT to the quantum do- 
main has recently attracted significant attention [7, 8]. 
Specifically, a quantum exchange FT, for the transfer 
of energy between two reservoirs maintained at differ- 
ent temperatures, has been derived in Refs. [4, 9, 10] 
using projective measurements, and in Refs. [11, 12], 
based on the unraveling of the quantum master equation 
(QME). These derivations assume that the interaction 
between the two thermal baths is weak, and can be ne- 
glected with respect to overall energy changes. Using the 
Keldysh approach, an exact analysis was carried out in 
[13]. However, it is valid only for harmonic systems. It 
is thus an open question whether a heat exchange FT is 
obeyed by an anharmonic quantum system strongly cou- 
pled to multiple reservoirs. 

From a practical point of view, understanding and 
controlling energy transport and heat dissipation in 
nanoscale junctions is crucial for making further progress 
in device miniaturization [14]. Theoretical studies adopt- 
ing simple models can reveal the role of different system 
parameters on the transport mechanisms [15-18]. How- 
ever, such treatments either assume weak coupling be- 
tween the nanoscale object and the environment, an as- 
sumption that is not always justified, or are limited to 
very simple models. 



It is our objective here to investigate quantum heat 
exchange in two-terminal impurity models: (i) To derive 
the SSFT for heat currents in open quantum systems, in- 
corporating anharmonic interactions, allowing for strong 
system-bath interactions ("strong coupling"), (ii) To ob- 
tain the CGF and gain explicit expressions for the heat 
current and its second moment, useful for understand- 
ing heat current characteristics for anharmonic-strongly 
coupled systems, (iii) To understand the role of non- 
markovian (memory) effects on the onset of the SSFT. 

Our analysis begins with a general model for the impu- 
rity, reservoirs and the interaction form. Describing the 
dynamics at the level of the noninteracting-blip approx- 
imation (NIB A) [19], a scheme accommodating strong 
system-bath interactions, we derive a QME for the sys- 
tem dynamics, under the markovian limit. Unraveling 
these equations into trajectories with a particular amount 
of net energy dissipated, e.g., to the R reservoir, a heat 
exchange SSFT is verified. We also obtain the CGF, 
independent of the particular physical realization. The 
scheme is exemplified on the two-terminal spin-boson 
model. In the nonmarkovian case a general symmetry re- 
lation is recovered, whereas the universal SSFT is reached 
in the markovian limit only. 

Model. — Consider a quantum impurity (system) 
placed between two thermal reservoirs (baths). No as- 
sumptions are made on the energy structure of the im- 
purity, thus anharmonic systems, with finite and uneven 
energy spacings, are comprised. Further, system-bath 
interactions are potentially strong relative to the system 
energetics. We adopt the dressed-tunneling Hamiltonian, 



H = J2^n\n)(n\+J2H, 

n v 

+ E (HHe-^ + \m){n\e in ™) , (1) 

n > rn 

where |n) denotes the impurity quantum states, cou- 
pled through the tunneling elements A nm , dressed by 
the baths operator fl nm = fl n mL + ^nmR- The opera- 
tors r2„ m „ depend on the coordinates of the v = L, R 
bath and may represent, for example, a collection of dis- 
placements or momentum operators as in the standard 
small polaron model [19]. Furthermore, different bath 
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operators may couple to different transitions. The ther- 
mal reservoirs H v are assumed to be in a canonical state, 
maintained at a temperature T v = fi^ 1 . Besides that, we 
do not specify the reservoirs, and they may be composed 
of fermions, spins, photons or phonons. The Hamiltonian 
(1) allows only for energy transfer processes between the 
two baths, mediated by a system excitation. Transfer of 
particles is not considered in the present study. 

population dynamics. — System dynamics is explored 
at the level of the NIB A scheme [20-22]: Applying the 
Born approximation [19] to the dressed Hamiltonian (1), 
equations of motion for the impurity reduced density ma- 
trix can be readily obtained [16]. This approximation is 
generally valid for A < lo c , where lo c is a cutoff of the 
reservoirs modes, at high temperatures and in the strong 
coupling regime [19]. Neglecting coherences and for sim- 
plicity, further applying the Markov approximation, we 
get quantum kinetic equations for the population p n , 



Pn = ~Pn ^ 



The transition rate from state n to m, C nrn (cu nrn ), is a 
convolution of L-induced and R induced processes [16], 

/oo 
e t0Jnm C nrn L(t)C nm R(t)dt 
-oo 

C n mL(^nm — ^)C nm R(uj)duJ. (3) 



-L 



Here LO nm = e n — e rn . The indices of C nm are ordered such 
that n > m. The i^-bath correlation function is given by 
the thermal average 



Cnmz/ (^) 



,(t) e -ifi„ 



.(0) 



(4) 



The operators are written in the interaction representa- 
tion, Q n mv(t) — e lHvt £l nmv e~' lH ' jt . In frequency domain 
we write C nmv (uj) = J_ dte lu)t C nmv (t), which are the 
elements in (3). As a result of microreversibility, detailed 
balance is satisfied for each reservoir, separately, 



G n rav{f-0) 
C nmv (—Lo) 



= e 



(5) 



Such a detailed balance condition does not hold for the 
combined rate C nm (to nm ) since it encloses both temper- 
atures through the bath-specific correlations C nmu . 

While the dynamics is simply described by a QME, it 
still encloses complex physical processes. Eq. (3) draws 
nontrivial transfer rates. For example, when the sys- 
tem decays making a transition from state n to m, it 
disposes the energy LO nm into both reservoirs coopera- 
tively; an energy to is dissipated into the R bath while 
the L bath gains (or contributes) the rest, Lo nm — lo. Sim- 
ilarly, excitation of the system occurs through an L-R 
compound process. We highlight the three non trivial 
mechanisms involved here, arising due to the strong cou- 
pling limit: (i) Non-resonance energy transfer processes 



are allowed, where each reservoir donates (absorbs) an 
energy which does not overlap with the system's energy 
spacings. (ii) Anharmonic processes are allowed. For ex- 
ample, in the context of vibrational energy transfer multi- 
phonon processes are incorporated within the relaxation 
rates C nm , see e.g., Eq. (14). (hi) The transport pro- 
cess takes place conjoining the reservoirs' dynamics in a 
non-additive manner, as discussed above. In contrast, 
the weak coupling limit, studied in [11, 12, 16] in the 
context of bosonic transfer, admits only resonant trans- 
mission processes and single phonon effects. Moreover, 
in the weak coupling limit the reservoirs additively act 
on the system [23]. 

Cumulant Generating function. — We define the func- 
tion P t (n,uj) as the probability distribution that within 
the time t a net energy lo has been dissipated into the R 
bath, with the system populating the n state at time t. 
For later use we also construct Pt{oo) = ^ n Pt{n, lo), the 
distribution of lo at t, irrespective of the system state. 
The time evolution of P t (n,uj) obeys 

/oo 
P t (m,Lo)C nmR (u - lo) 

xC nmL {ui - LO - LO nm )du} 

/oo 
C nm R(u)C nmL (uj nm - u})dui.(6) 

This can be justified by energy-resolving the popula- 
tion dynamics in (2), then collecting the matching en- 
ergy terms from the left side and the right side of the 
equation. The first term here describes a process where 
by the time t a net energy ui has been damped into R, 
whereas the system occupies the state m. At the moment 
t the system (assisted by the bath) transits from m — > n, 
further dissipating an energy lo — uj into the R reser- 
voir. Similarly, the second term collects all transitions 
which deplete P t (n,Lo). Next we introduce the counting 
field x an d Fourier transform the resolved probabilities, 
Pt(n,x) = J^dLoe^xPtin,^), yielding 

Pt(n,X) = -Pt(n,X) ^2 C nm(Unm) 

+ ]T Pt(m,x)/+ n (x)+ E Pt(™,x)fnm(x)- (7) 

m>n m<n 

For brevity, we introduce the short notation 

/oo 
e^C nmR (Lo)C nmL (±uj nm - Lo)dLO. (8) 
-oo 

These equations can be encapsulated in a matrix form 
l*(X)*)) = -Kx)\^(x,t))i with * a vector of the prob- 
abilities P t (n,x)- We define the characteristic func- 
tion Z(x,t) = (I\^(x,t)), with (J|, as a left vector of 
unity, and the cumulant generating function G(x) — 
lhm-j.oo j\nZ(x,t), recovered as the negative of the 
smallest eigenvalue of the matrix fi. 
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Steady-state fluctuation theorem. — We now prove that 
G(x) = G(iA(3 - x), implying that a SSFT for heat ex- 
change holds. In order to derive this relation we analyze 
the symmetry properties of the matrix fi. For clarity, we 
explicitly write it for a three-state impurity 

/ -ftiix) -/ 3 + i(x)\ 

Kx) = -/ 2 "i(x) ^2,2 -/ 3 + 2 (x) (9) 

V -/ 3 1(X) -/ 32 (X) M3,3 / 

The diagonal terms /z^, constitute the decay rates from 
each level, and are independent of \. The characteristic 
polynomial D^( X )(A), with the roots A, is given by 

D Kx) (X) = / 3 + !(x) [/ 2 "i(x)/32(x) - (A - M2, 2 )/ 3 - 1 (x)] 

-/ 2 + i(x) [/ 2 "i(x)(A - M 3i3 ) - / 3 + 2 (x)/ 3 l(x)] 

+(A - Aii.i) [(A - M2;2 )(A - M3 , 3 ) - / 3 - 2 (x)/ 3 + 2 (x)] • 

One can show that the following three properties hold: 
(i) fi(x) is symmetric under the operation fnmix) ~ ► 
fnmix)- Thus, the roots A are also symmetric in this 
respect, (ii) Each element in the characteristic polyno- 
mial is cyclic, in the sense that a series of transitions 
must end at the initial state. For example, the product 
/ 3 ~i(x)/2i(x)/3~2(x) describes a relaxation process from 
state 3 to 1, followed by an excitation from state 1 to 
2, finishing with an excitation term /^(x); bringing the 
system back to state 3. (iii) The correlation function 
fnmix) satisfies the identity 

/+ ro (iA/3- X ) = e^-/- ro (x), (10) 

gathered by manipulating Eq. (8) with (5). Under these 
three properties we prove that D^i x \{X) = D^ i /^p_ x ^{X): 
The symmetric terms in the characteristic polynomial are 
mapped one onto the other as a result of the symme- 
try (10) whereas the system dependent prefactors, i.e., 
the term e^ 1 -"™"* in Eq. (10) overall cancel, a result of 
the cyclic property (ii). We conclude that the eigenval- 
ues of ft satisfy a symmetry relation, and in particular 
G(x) = G(iA(3 — x). The probability distribution of 
co is obtained as Pt(oj) = j- dxZ(x,t)e~ lxu> . Since 

Z( X ,t) ~ e G ^* in the long time limit, a heat exchange 
fluctuation relation is resolved 

i im im4M.^. (11) 

We emphasize: This relation has been derived without 
specifying neither the system energy structure and its 
interaction with the reservoirs, nor the form of the reser- 
voirs. It allows for strong coupling between the impu- 
rity and the baths, reflected in the transition rates C nm , 
mixing L-R processes in a non-additive manner. More- 
over, an explicit expression for the CGF, G(x), can now 
be written, bearing analytical expressions for the current 
cumulants, as we achieve below for the spin-boson model. 

Spin boson model. — The equilibrium spin-boson (SB) 
model, referring to a spin immersed in an equilibrated 



boson reservoir, is an eminent model in chemistry and 
physics, useful for describing, e.g., solvent assisted elec- 
tron transfer reactions and the Kondo resonance [19]. 
The nonequilibrium spin-boson model, where the spin 
is coupled to more than one thermal reservoir, has been 
suggested as a prototype model for exploring heat trans- 
fer through nanojunctions [16, 18]. We now analytically 
obtain the CGF, thus the current and its moments, for 
the nonequilibrium SB model at strong coupling, 

H = Y^ + ^CTx+CT z ^Aj, I/ (6] )1/ + 6j>)+^a;j6t i/ 6 J >. 

(12) 

Here a x and a z are the Pauli matrices, ojq is the energy 
gap between the spin levels, and A is the tunneling en- 
ergy. The two reservoirs include a collection of uncoupled 
harmonic oscillators, bj v (bj tV ) is the bosonic creation 
(annihilation) operator of the mode j in the v reservoir. 
The parameter accounts for the system-bath inter- 
action strength. The Hamiltonian is transformed to the 
displaced bath-oscillators basis using the small polaron 
transformation [19], H s = U^HU, U = e ia * n / 2 , 

H S = fa z + | (a + e in + a.e~ in ) + £ Uj b] >v b j>v (W) 

where a± = \ {p x ^io~y) are the auxiliary Pauli matrices, 
= 52 v Sl v , and n v = 2iJ2j - &,>)• Under the 

NIBA, the system population obeys a convolution-type 
master equation [20-22] ((cr z ) = p\ —po), 

Pi = / e- Q ' (t - s) cos[ Wo (t - s) - Q"(t - s)} Pl ( S )d S 

* Jo 

+ 4p / e- Q '(*- s ) cos[wo(f - s) + Q"(t - s)}p (s)ds, 

* Jo 

with conserved total occupation po(t) + pi(t) = 1. The 
function Q(t) = ^ v Q v (t), made of a real and imaginary 
components, Q v (t) = Q' v {t) +iQ"(t), is defined by 

Q'„(t) = [™ l^[l-cos(ut)][l + 2n v (u)]du, 
Jo koj 2 

Q'l(t) = ^ ^-sm(ojt)dw. (14) 
Jo ^ 

Here J v (uj) = Air ^ . )^ v 8{oj — uoj) is the j/-bath spec- 
tral function, n v (oj) is the Bose- Einstein distribution. 
We have carried out the analysis in the nonmarkovian 
limit by generalizing Eq. (6), to describe the dynamics 
of Pt(n,u>L,WR), for the transfer of uj v net energy to the 
v bath by the time t. Introducing two counting fields 
Xi.2, then following the procedure outlined in Ref. [24] 
(applying Fourier transform and Laplace transform on 
the resolved equation of motion, analyzing the poles of 
the resolvent), we can prove that the CGF satisfies [25] 

G(x 1 ,x») = G(ip L -xi,0R-X»)- (15) 
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Only in the markovian limit the symmetry is given in 
terms of the affinity as G(x) = G(iA/3 — x). Thus, while 
microreversibility is sufficient for deriving the basic sym- 
metry relation (15), the SSFT holds only under more re- 
strictive conditions, dictated here by the bath relaxation 
timescale [7, 9]. 

In the markovian case the QME for the population dy- 
namics [Eqs. (2)-(3)] follows pi = — C(uJo)pi + C(— ujo)po: 
with the rates C(w ) = e iuiot C L {t)C R {t)dt- C v {t) = 
e -Qv(t) 4 Since only a single correlation function matters, 
the level indices were discarded. Following Eqs. (6)-(8), 
we identify the matrix (1 by 



A(x) 



C(-wo) -f + {x) 

-f-(x) c(w„) 



(16) 



with f ± (x)= J_^ o e i " x C fl (w)C* L (±w -a;)dw. Its small- 
est eigenvalue is 



G(x) = -~[CM+C(-a; )] 



+ 2 #W - C-(-^o)) 2 + 4/-(x)/+(x)- (17) 
The averaged heat current can be readily obtained, 




FIG. 1: Nonequilibrium spin-boson model: Plot of Pt(w) at 
various times. The inset demonstrates the validity of the 
SSFT. T L = 3, T R = 2, Er = 1, uo = 0.5, t = 20 (o), 
t = 100 (dotted) and t = 400 (□). 



(J) 



(u) t _ dG(x) 



t d(i X ) 
C r (-u>)C l (-ujo + w)p 
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/ — OO 



Cr(u)Cl{uo - w)pi 
(18) 



The population here is calculated in steady-state, p — 
C (luq) / [C (loq) + C(— u>o)]- This expression was heuristi- 
cally suggested in Rcf. [16], and here it is derived from 
the basic dynamics. Note that the details of the function 
Q(t) are not utilized in this derivation. Furthermore, the 
averaged current stays intact for nonmarkovian systems 
[24] . The formal structure for the noise power is given by 



(S) 



d 2 G(x) 



d{i X f 



= -2 



C(w ) + C(-wo) 



ojC-(uj)duj I ujCj r {uj)duj + (jy 



dujui^ 



(19) 



where we defined C±(w) = C_r(±w)Cl(±o;oT w )- We can 
also plot the distribution Pt(oo). Assuming high temper- 
atures T„ > wo and strong coupling, Eq. (14) can be 
simplified, Q' v {t) = E»TJ 2 , Q'l(t) = E v T t, with the re- 
organization energy defined as — 4\? v /u>j [26]. 
Using this form, Fig. 1 displays the entropy production 
distribution and the validity of the SSFT (inset). 

To conclude, a heat exchange SSFT has been derived 
for quantum systems incorporating strong system-bath 
interactions and anharmonic effects. Our study provides 
closed expressions for the CGF, useful for deriving the 
distribution of heat fluctuations, the averaged current 
and the thermal noise power. For the spin-boson model 
one can show that in the nonmarkovian case the SSFT 
does not generally hold. It is satisfied in the marko- 
vian limit, when energy conservation is enforced. Future 
work will be devoted to generalizing our study to systems 
showing coherence effects. 
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